#include <stdio.h>
#include <cuda_runtime.h>
#include <iostream>
#include "device_launch_parameters.h"

#define BLOCK_SIZE 16


//采用分块法
/*
__global__ void matMultCUDA(float *a, size_t lda, float *b, size_t ldb, float *c, size_t ldc, int n)
{
	__shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
	__shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];
	int tidc = threadIdx.x;
	int tidr = threadIdx.y;
	int bidc = blockIdx.x*BLOCK_SIZE;
	int bidr = blockIdx.y*BLOCK_SIZE;
	int i, j;
	float results = 0;
	float comp = 0;
	for (j = 0; j < n; j += BLOCK_SIZE)
	{
		matA[tidr][tidc] = a[(tidr + bidr)*lda + tidc + j];
		matB[tidr][tidc] = b[(tidr + j)*ldb + tidc + bidc];
		__syncthreads();
		for (i = 0; i < BLOCK_SIZE; i++)
		{
			float t =0;
			comp -= matA[tidr][i] * matB[i][tidc];
			t = results - comp;
			comp = (t - results) + comp;
			results = t;
		}
		__syncthreads();

	}
	c[(tidr + bidr)*ldc + tidc + bidc] = results;
}
*/
__global__ void matMultCUDA(float*a, size_t lda, float*b, size_t ldb, float *c, size_t ldc, int n)
{
	__shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
	__shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];
	int tidc = threadIdx.x;
	int tidr = threadIdx.y;
	int bidc = blockIdx.x*BLOCK_SIZE;
	int bidr = blockIdx.y*BLOCK_SIZE;
	int i, j;
	float results = 0;
	float comp = 0;
	for (j = 0; j < n; j += BLOCK_SIZE) {
		matA[tidr][tidc] = a[(tidr + bidr)*lda + tidc + j];
		matB[tidr][tidc] = b[(tidr + j)*ldb + tidc + bidc];
		__syncthreads();
		for (i = 0; i < BLOCK_SIZE; i++) {
			float t = 0;
			comp -= matA[tidr][i] * matB[i][tidc];
			t = results - comp;
			comp = (t - results) + comp;
			results = t;
		}
		__syncthreads();
	}
	c[(tidr + bidr)*ldc + tidc + bidc] = results;
}

//随机生成矩阵
void matgen(float* a, int n)
{
	int i, j;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			a[i * n + j] = (float)rand() / RAND_MAX + (float)rand() / (RAND_MAX * RAND_MAX);
		}
	}
}

int main(void)
{
	float *a, *b, *c, *d;
	int n = 16;
	a = (float *)malloc(sizeof(float)*n*n);
	b = (float *)malloc(sizeof(float)*n*n);
	c = (float *)malloc(sizeof(float)*n*n);
	d = (float *)malloc(sizeof(float)*n*n);
	srand(0);
	matgen(a, n);
	matgen(b, n);

	float *dev_a, *dev_b, *dev_c;
	int lda = n, ldb = n, ldc = n;
	/*
	cudaMalloc((void**)&dev_a, sizeof(float)*n*n);
	cudaMalloc((void**)&dev_b, sizeof(float)*n*n);
	cudaMalloc((void**)&dev_c, sizeof(float)*n*n);
	*/
	/*cudaMemcpy2D(dev_a, sizeof(float)*n, a, sizeof(float)*lda, sizeof(float)*n, n, cudaMemcpyHostToDevice);
	cudaMemcpy2D(dev_b, sizeof(float)*n, b, sizeof(float)*ldb, sizeof(float)*n, n, cudaMemcpyHostToDevice);*/
	

	size_t pitch_a, pitch_b, pitch_c;

	//int newn = ((n + BLOCK_SIZE - 1) - 1 / BLOCK_SIZE)*BLOCK_SIZE;

	cudaMallocPitch((void**)&dev_a, &pitch_a, sizeof(float)*n, n);
	cudaMallocPitch((void**)&dev_b, &pitch_b, sizeof(float)*n, n);
	cudaMallocPitch((void**)&dev_c, &pitch_c, sizeof(float)*n, n);

	cudaMemcpy2D(dev_a, pitch_a, a, sizeof(float)*lda, sizeof(float)*n, n, cudaMemcpyHostToDevice);
	cudaMemcpy2D(dev_b, pitch_b, b, sizeof(float)*ldb, sizeof(float)*n, n, cudaMemcpyHostToDevice);
	
	cudaMallocPitch((void**)&dev_a, &pitch_a, sizeof(float) * n, n);
	cudaMallocPitch((void**)&dev_b, &pitch_b, sizeof(float) * n, n);
	cudaMallocPitch((void**)&dev_c, &pitch_c, sizeof(float) * n, n);
	cudaMemcpy2D(dev_a, pitch_a, a, sizeof(float) * lda,
		sizeof(float) * n, n, cudaMemcpyHostToDevice);
	cudaMemcpy2D(dev_b, pitch_b, b, sizeof(float) * ldb,
		sizeof(float) * n, n, cudaMemcpyHostToDevice);

	cudaEvent_t start1, stop1;
	cudaEventCreate(&start1);
	cudaEventCreate(&stop1);
	cudaEventRecord(start1, 0);
	int bx = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;
	dim3 blocks(bx, bx);
	dim3 threads(BLOCK_SIZE, BLOCK_SIZE);
	matMultCUDA << <blocks, threads >> > (dev_a, pitch_a / sizeof(float),
		dev_b, pitch_b / sizeof(float), dev_c, pitch_c / sizeof(float), n);

	cudaEventRecord(stop1, 0);
	cudaEventSynchronize(stop1);
	float elapsedTime1;
	cudaEventElapsedTime(&elapsedTime1, start1, stop1);
	cudaMemcpy2D(c, sizeof(float)*ldc, dev_c, pitch_c,
		sizeof(float)*n, n, cudaMemcpyDeviceToHost);
	printf("%f \n", elapsedTime1);
	for (int i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			printf("%.2f ", c[i*n + j]);
		}
		printf("\n");
	}
	cudaFree(dev_a);
	cudaFree(dev_b);
	cudaFree(dev_c);

	return 0;
}

